function [pr,CS]=pr_fun(beta,X,TC_idx)
KX=size(X,2); L=length(TC_idx);
% Build vectors of expressions containing parameters
delta=X*beta;  ed=exp(delta); % exponential of delta
if sum(isinf(ed)) >0
   fprintf('pause since delta is inf');
    pause
end
sed=my_accumarray(TC_idx,ed); %sum of ed(j) within group TC
pr=ed./(1+sed); 
if nargout==1, return; end
CS=log(1+sed)./beta(1);
end
